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Abstract 

Evaluation of pfaffians arises in a number of physics applications, and for some 
of them a direct method is preferable to using the determinantal formula. We 
discuss two methods for the numerical evaluation of pfaffians. The first is tridi- 
agonalization based on Householder transformations. The main advantage of 
this method is its numerical stability that makes unnecessary the implementa- 
tion of a pivoting strategy. The second method considered is based on Aitken's 
block diagonalization formula. It yields to a kind of LU (similar to Cholesky's 
factorization) decomposition (under congruence) of arbitrary skew-symmetric 
matrices that is well suited both for the numeric and symbolic evaluations of 
the pfafhan. Fortran subroutines (FORTRAN 77 and 90) implementing both 
methods are given. We also provide simple implementations in Python and 
Mathematica for purpose of testing, or for exploratory studies of methods that 
make use of pfaffians. 
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1. Introduction 

In a number of fields in physics, the formal equations derived from the theory 
make use of the pfaffian of some skew-symmetric matrix appearing in the theory. 
For example, the pfaffian arises in the treatment of electronic structure with 
quantum Monte Carlo methods [1], the description of two-dimensional Ising spin 
glasses [2] J s-iid the evaluation of entropy and its relation to entanglement 
Pfaffians occur naturally in field theory and nuclear physics in formalisms based 
on fermionic coherent states 3, [1,0, 01 • A recent application is to the overlap 
of two Hartree Fock Bogoliubov (HFB) product wave functions [^], needed for 
nuclear structure theory. While there is a simple formula for the pfaffian of a 
skew-symmetric matrix M in terms of the determinant. 



the so-called "sign problem of the overlap" [9] associated with the square root 
motivates the use of numerical algorithms that evaluate it directly. The most 
straightforward method, the rule of "expanding in minors" [10], has bad scaling 
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with the size of the matrix and is prohibitive for large matrices. In this paper 
we discuss two alternative methods that have the same scaling property as the 
normal algorithms for the determinant. The methods are implemented in 
the FORTRAN 77 and 90 subroutines provided in the accompanying program 
library. We also comment on the practical implementation of the two methods 
in Mathematica and in the Python programming language. 



2. Evaluation of the PfafRan 



The Pfaffian pi{A) is reduced to a simple form that is easily evaluated by 
making repeated use of transformation formula given in [Appendix A[ 

pi{B'AB) = det(B)pf(A). (2) 

In order to perform the numerical evaluation of the Pfaffian of a complex skew- 
symmetric matrix A we reduce the skew-symmetric matrix to a tridiagonal form 
Atr by using unitary matrices U. Once it is in this form, the evaluation of the 
pfaffian is straightforward (see below). 



2.1. Reduction to tridiagonal form by mean of Householder transformations 

In this method, we will use the well-known Householder transformations 
(m to reduce A to tridiagonal form. We present it in some detail because the 
generalization to the complex number field is not entirely trivial. 

Complex Householder transformations have the form 



1-2 



u(E) u~' 



(3) 



where u is an arbitrary complex complex row vector u — (ui, M2, . . . , un) and 
(u (8) ~ UiU*. The vector u must be chosen to zero all the elements of a 

vector X except a given one. If we take u = e^'^'^^^^^lxjej, with (ej)^ — Sjk, 
it can be easily proved that 



as required. The freedom on the sign in the expression defining the vector u can 
be used to make sure that the vector u is non zero. The rest of the Householder 
tridiagonalization procedure follows exactly as in the real case. Consider a 
skew-symmetric matrix of dimension N (even) 
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-ai2 






-ai3 
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The Householder transformation matrix is 



/ 1 



Pi = 



V 



J 



where -'^'Pi is buih by using Eq. ^ and taking the vector x (of dimension 
N-1) as (ai2,ai3, ■ • ■)'^ ■ The resulting transformed matrix is given by 



/ 



V 






ki 


-fcl 










where fci = ±e' ^"^8(012) ^j^^ fj^^ matrix -^^A is skew-symmetric and given 
by = iN'i)p^{N-i)j!^(N-i)pT^ Performing this procedure a total of N-2 

times we end up with a tridiagonal and skew-symmetric matrix 



/ 





p 



N-2 ■ 



..P^PiAP^P^ 



fcl 







fc2 



\ 






-fcjv- 



fcw- 




V : ■ -kN-i J 

(5) 

Using now a known property the Pfaffian (see [Appendix A| we can deduce from 
the above identity that det(Pi) . . . det(P/v_2)pf(^) = pf(^Tfl) where Atr is the 
triagonal and skew-symmetric matrix of the right hand side of Eq. ^ . Taking 
into account that the determinant of any Householder matrix is -1 and that N is 
even, we can express the Pfaffian of A in terms of the pfafhan of the tridiagonal 
Atr 

pf(A) = pi{ATR) 

As will be shown below the Pfaffian of a tridiagonal skew-symmetric matrix is 
simply given by fcifcs . . . fcAr-i = 77 ^]^ fc 2z-i (this result can also be obtained 
using the "minor expansion" formula |lO| ) and finally we obtain 



pf(A) 



N/2 

II '^2 

1=1 



(6) 



In terms of numerical stability, the Householder transformation is very ro- 
bust and there is no need to consider any "pivoting" strategy common to other 
methods. However, the presence of the square root of x and the argument 
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arg(xj) of complex quantities prevents an easy implementation of the House- 
holder tridiagonalization procedure for symbolic computation. For this purpose 
the second method described in the next section is far easier to implement. 

2.2. Aitken's block diagonalization formula 

There is an alternative method for the calculation of the pfaffian, which is 
also well suited for a symbolic implementation and that relies on an expression 
for the pfaffian of a bipartite skew-symmetric matrix. Let us start with a general 
skew-symmetric matrix A (dimension N, even) given by 

where R and S are square skew-symmetric matrices and Q is a general rectan- 
gular matrix (to account for the case where R and S have different dimensions). 
Using Aitken's block diagonalization formula (see for an early use of the 
formula and [l3| for a recent and thorough reference) for a bipartite matrix we 
obtain 

I Q \ ( R Q \ ( I -R-^Q 
Q^R-' I ) [ -Q^ S ) [ I 

R 



S + Q^R-^Q 



(8) 



where the matrix S + Q^R~^Q is referred to in the literature as the Schur 
complement of the matrix A (see, for instance, [l3|)- For the special case of a 
skew-symmetric matrix A, the matrices R and S are also skew-symmetric and 
the transformation of the matrix A is a congruence (i.e. the matrix acting on 
the left hand side of A is the transpose of the one acting on the right hand side). 
Denoting 

P^=( nri-i M (9) 



Q'R-' I 



Eq. (ID) becomes 



^i^^i - 1^ S + Q^R-'Q 
An equivalent expression involving S^^ instead of R^^ is easily obtained 



P^AP^ = 



R + QS'^Q^ 
S 



with 



C ) (-) 

By using the property pilP'^AP) — det(P)pf(A) (see [Appendix A| and taking 
into account that detPi — detP2=l, we come to 

pf(^) = pf(P)pf(5 + Q^p-iQ) (11) 
= pf(i? + Q5-iQ^)pf(5) (12) 
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Another nice property of the matrices Pi and P2 is that their inverses can be 
obtained very easily 



and 



' - { I ' ) (14) 



These expressions of the inverses explicitly show that both Pi and P2 are not 
orthogonal matrices. 

Let us now apply the above result to an arbitrary skew-symmetric matrix of 
dimension N — 2M which is written in block form as 

/ An-i An \ 

A = -Al_^ ajv-i,jv (15) 

V -O^N-lM / 

where A'^^^ is a skew-symmetric square matrix of dimension N — 2 = 2{M — 1) 
and Aff-i and Aat are column vectors An-i — {A^jv-i, i = 1,A^ — 2} and 
An — {Ai N, i — 1, N — 2} both of dimension {N — 2) x 1. In the language of 
Eq ([7|) the matrix R is the matrix A^^\ the matrix Q is a rectangular matrix 
of dimension 2 x (N — 2) made of the two column vectors, A^-i and An and 
finally the matrix S is the 2x2 skew-symmetric matrix with matrix element 
iS'12 = ciN-i,N- Using the ideas of Aitken's block diagonalization formula, it is 
easy to shows that the matrix A = DfADi is in block diagonal form 

A^'^ 

aN-i,N I (16) 

-aN-i,N 

with a matrix Di of the form 

Di^\ X 10 1 (17) 



where lAr_2 stands for the identity matrix of dimension N — 2 and both X and 

Y are row vectors of dimension 1 x [N — 2) and given by AT = — a^"'^-^ and 

Y = a~^_^ j^A^_^. In the above equation 1161 the skew-symmetric matrix A*-^' 
is given by 

i(i) = ^(1) + AN{aN-i,Nr^Al_i - AN-ia-^\NA^N (18) 

Taking into account that detZ^i — 1 then pf(A) — pf(A) = a^r^ ijvpf(A*^^-'). 
The algorithm can be applied recursively to A^^^ to obtain 

pf(j4) — aN-i.N'i'N-z Ar-9Pf(^*'^'') so as to reduce, after M — 1 iterations, the 
computation of the pfafhan to the product of the corresponding elements. 
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This procedure can be easily implemented for a skew-symmetric tridiagonal 
matrix, as the transformed matrices in Eq ([T8| coincide with the original ones; 
for instance, A^^^ = A^^\ As a consequence, the pfafBan of a tridiagonal matrix 
is given by 

/ 



pf 






di 


















d2 













-d2 






































d2N 













— d2N-l 






N 



dids . ..d2N-i = 



i=l 



2.2.1. Pivoting 

As a consequence of the division by matrix elements like aN-i,N in the first 
iteration, the numerical stability of the algorithm requires the use of pivoting 
strategy in the implementation of the method. Full pivoting amounts to search 
the whole matrix for the matrix element with the largest modulus and exchange 
it with the required matrix element. For instance, in the first iteration of the 
procedure, the matrix element ap^q {p < q) with the largest modulus is searched 
for and exchanged with the matrix element aN-i,N- In this way we avoid 
dangerous divisions by small (or even zero) matrix elements. We have to take 
into account that in the present case, the exchange of both columns and rows 
is required to preserve the skew-symmetric nature of the matrices involved. To 
carry out the exchange of rows and columns we will use the exchange matrix 
P{ij) that, when applied to the right of an arbitrary matrix, exchanges columns 
i and j. The exchange matrix is given by the matrix elements 

P{ij)ki = Ski - Si^iSi^k - 5j,iSj,k + Si^iSj^k + S],iSi,k. (19) 

To exchange the corresponding rows we have to apply P{ij)'^ to the left of the 
matrix (notice that P{ij) is symmetric). With the help of these matrices we 
can write the matrix after pivoting ap^q with aN-i.N (and aq^p with aN,N-i) as 

Ap = P'^iN - l,p)P^(iV, q) A P{N - l,p)P{N, q) 

As a consequence of such exchange and taking into account that det P(ij) = — 1 
we can conclude that the pfafhan of A does not change by the pivoting procedure. 
Finally we obtain 

A = P{N, q)P{N - l,p)ApP^{N - l,p)P'^{N, q) 

= PiN, q)PiN - l,p)D'[-'ApD-'P^{N - l,p)P^(iV, q) 

where Ap has the same structure as S in Eq. (I16p . As before, pi{A) — pi{Ap) = 
{Ap)j^_i pf(^p'') and repeating recursively the whole procedure M — 1 times 
we obtain the pfafhan as the product of the corresponding matrix elements. 
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2.2.2. Cholesky like decomposition of a skew- symmetric matrix 

Although it is not necessary in order to compute the pfafBan, it can be useful 
to show that even with pivoting we can write the matrix A as 

A = PL^ALP (20) 

where P is the product of exchange matrices as in Eq p^ . L is the product of 
matrices of the D^^ type, Eq (fT7|) . and therefore is a lower triangular matrix 
with ones in the main diagonal and finally, A is a skew-symmetric matrix in 
canonical form, i.e. a block diagonal matrix with skew-symmetric, 2x2 blocks 
in the diagonal. This decomposition of a general skew-symmetric matrix A 
resembles the Cholesky decomposition of a general matrix and can be useful in 
formal manipulations like, for instance, the inversion of the matrix A. In order 
to show that Eq (|20p holds the only required property is that, when applying the 
pivoting procedure to A^p^ the exchange matrices required P{N—2, s)P{N—i, r) 
have the property of preserving the structure of the matrix Z?i(and its inverse). 
For instance, 

dI -^P{N - 2, s)P{N - 3, r) = P{N - 2, s)P{N - 3, r)Dj 

with ~^ a matrix that is obtained from ~^ by exchanging rows N — 2 and 
s and rows — 3 and r and therefore has the same upper triangular structure 
with ones in the diagonal as the original matrix dJ . Using this property 
we can move all the exchange matrices to the right (or to the left) and the 
remaining matrix will be the product of triangular matrices (lower for products 
involving D~^) with ones in the diagonal. 

As mentioned earlier, Aitken's method is better adapted to symbolic eval- 
uations. However, one must take care that in each step of the process some 
specific matrix elements are non-zero. 

3. Fortran implementation 

The implementation of the algorithms considered in this paper in a high 
level computer language is straightforward. However, specific code in FOR- 
TRAN (both 77 and 90, real and complex arithmetic) is provided along with 
this paper. The algorithms are easy to follow and the comments included in 
the code are useful guides. Just a few comments are in order: to implement 
the tridiagonalization procedure in Fortran, it is advantageous to use the BLAS 
package [14i] to perform the required matrix by vector multiplication and rank 
two update. Unfortunately there are no equivalent in the skew-symmetric case 
of the routines SYM (to multiply a symmetric matrix by a vector) or SYR2 
(to perform a symmetric rank two update) but the general procedures GEMV 
and GERU can be used instead. On output, both the pfaffian of the matrix 
and the set of vectors required to bring it to tridiagonal form are returned. 
For the implementation of the method based on Aitken's block diagonalization 
formula a pivoting strategy is required. We have used full pivoting in our imple- 
mentation due to its robustness. The routines provided only require the upper 
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part of the skew-symmetric matrix. The lower part is destroyed and replaced 
with the tridiagonal transformation matrix that brings the skew-symmetric ma- 
trix to canonical form upon congruence. An integer vector is also returned to 
reconstruct the required exchange of rows and columns. 

Perhaps the best test to check the validity of the two implementations is 
to compute the pfaffian of a skew-symmetric matrix using both procedures in 
order to compare the output. If it is the same up to a given accuracy then 
it is very likely that the two implementations are correct. We have writen a 
test program (also included in the distribution) that generates skew-symmetric 
matrices of given dimension with random entries and compute the pfafRan using 
both techniques. In our tests the pfaffians computed both ways coincide up to 
one part in 10^" with dimensions of the matrices of one thousand. This result 
also supports the adequacy of the implementation in terms of numerical stability. 
Another possibility to test the numerical implementation is to use the analytical 
formula given in |Appendix B| for a specific kind of 8 x 8 matrices. A test program 
implementing this approach has also been included in the distribution. 

To finish this section we will briefly comment on the timing of the FOR- 
TRAN numerical implementations mentioned. In a modern personal computer 
under Linux the computation of the pfaffian of a 100 x 100 matrix takes a few 
milliseconds in both implementations and the timing scales roughly as the cube 
of the dimension of the matrix in such a way that for matrices of a 1000 x 1000 
dimension the time is of the order of a few seconds. 

4. A simple Python implementation 

We provide here a simple implementation of the tridiagonal reduction method 
(see [12] and [15]) in Python, which may be useful for testing purposes. It is 
similar to the Householder, but it only use simple row and column operations 
that have determinants of unity. The code is: 

from numpy import * 
def pfaff_py(m) : 

mat=copy (m) 

ndim = shape (mat) [0] 

tl=1.0 

for j in range (ndim/2) : 
tl *= mat [0,1] 
print 'tl', tl 
if j <ndim/2-l : 

ndimr=shape (mat) [0] 
for i in range(2,ndimr) : 
if mat [0,1] != 0.0 : 

tv=mat [1 , : ] *mat [i , 0] /mat [1,0] 
mat [i , : ] -= tv 
tv=mat [ : , 1] *mat [0 , i] /mat [0 , 1] 
mat [ : , i ] -= tv 
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else : 

print 'need to pivot' 
raise Exception 
mat=mat [2: ,2:] 
return tl 

The user should be cautioned that the algorithm is not guaranteed to be stable 
without an additional pivot step. Also, the matrix is assumed to have been 
constructed with the array function in the Numpy library. 

5. A simple Mathematica implementation 

We also provide a simple Mathematica implementation of the method based 
on Aitken's block diagonalization formula. As mentioned above, this method 
requires pivoting to avoid divisions by small (or zero) numbers. In the symbolic 
implementation, this issue is solved by replacing the denominator by a variable 
(00 on the implementation below) in case it is zero and an additional limit 
when the variable tends to zero is performed at the end. The two Mathematica 
modules required are: 

Aitken [M_ , n_ , 00_] : = 
Module [{MM=M,i,j,p>, 

If [MM[[n-l,n]]==0,MM[[n-l,n]]=00;MM[[n,n-l]]=-00] ; 
p=MM[[n-l,n]] ; 
For [i=l , i<=n-2 , i++ , 
For [j=l,j<=n-2,j++, 

MM[[i,j]]=M[[i,j]]+(MM[[i,n]]*MM[[j,n-l]]-MM[[j,n]]*MM[[i,n-l]] 

] 

]; 

MM] ; 

pfaffian[S_] := 
Module [{T=S,n,p}, 
n=Length [T] /2 ; 

If [T [ [2*n- 1 , 2*n] ] ==0 , T [ [2*n- 1 , 2*n] ] =00 ; T [ [2*n , 2*n- 1] ] =-00] ; 
For [n=Length [T] /2 ; p=T [ [2*n- 1 , 2*n] ] , n> 1 , n- - , 

T=Aitken[T,2*n,00] ; 

p=p*T[[2*(n-l)-l,2*(n-l)]] 

]; 

Limit [p,00->0]] ; 

6. Conclusions 

The issue of how to compute both numerically and symbolically the pfaflian 
of a skew-symmetric matrix has been addressed using two different approaches. 
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Numerical stability issues are discussed and methods to assure the desired ac- 
curacy are fully incorporated. A collection of subroutines and test programs 
in FORTRAN (both 77 and 90, double precision and complex) are provided. 
A few comments on the implementation of the algorithms in Mathematica and 
Python are also given. 
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Appendix A. Definition and basic properties of the pfafflan 

The pfaSian of a skew-symmetric matrix R of dimension 2^^^ and with matrix 
elements is defined as 

Pf(^) = e(-P)^'^l^2»'^3^4^^5^6 • ■ • ?'2n-l,2n 

Perm 

where the sum extends to all possible permutations of ii, . . . , i2n and e(P) is 
the parity of the permutation. For matrices of odd dimension the pfaffian is 
by definition equal to zero. As an example, the pfaffian of a 2 x 2 matrix R 
is pf(i?) = ri2 and for a 4 x 4 one p{{R) = ri2'f34 — T'i3r24 + ri4r23- Useful 
properties of the pfaffian are 

p{{P'^RP) = det(P)pf(ii), (A.l) 

Pf( -R^ )=(-l)^(^-^'/^det(i?) 
where the matrix R is N x N and 

Pf ( ^2 ) = Pf(^l)Pf(^2) 

where Ri and i?2 are skew-symmetric matrices. The matrices may be defined 
on the real or on the complex number fields. 

Appendix B. Pfaffian of a test matrix 

In this appendix we give the expression of the pfaffian of a test matrix which 
is big enough as not to be trivial but on the other hand is small enough as to 
render the explicit expression of the pfaffian manageable. The expression given 
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below can be used to check both numerical and symbolic implementations of 
the pfaffian. 

Consider the two general skew-symmetric matrices of dimension 4 



M = 



and 



N 



1 





fi 


mil 


^7112 \ 




-/i 





m2i 


m22 




-TOii 


-m2i 
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-mi2 


-m22 


-/2 
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/ 





51 


nil 


ni2 \ 




-.91 





^^21 


?722 




-nil 







52 
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-ni2 


-n22 


-52 


/ 



where the matrix elements can be complex numbers. With these two matrices 
and the identity 4x4 matrix we build the skew-symmetric matrix 



S = 



N -I 
I -M* 



of dimension 8x8 (see Ref [8| for the physical context of this matrix). It is 
relatively easy to compute its pfaffian 

pf[5] = 1 + /i5i + /2*52 + rnliTiii + ■m*22n22 + m*i2ni2 + mlin2i + 

+ (/r /2* - "^iim22 + ml2mli){gig2 - niin22 + ni2n2i) 
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